Local sublattice-symmetry breaking in rotationally faulted multilayer graphene 
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Interlayer coupling in rotationally faulted graphene multilayers breaks the local sublattice- 
symmetry of the individual layers. We present a theory of this mechanism, which reduces to an 
effective Dirac model with space-dependent mass in an important limit. It thus makes a wealth of 
existing knowledge available for the study of rotationally faulted graphene multilayers. We demon- 
strate quantitative agreement between our theory and a recent experiment. 
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Experiments indicate that the 10-100 individual 
graphene layers grown on the carbon-terminated face of 
SiC arc surprisingly well decoupled from one another 
electronically. Early spectroscopic measurements 0, 0] 
found a linear low-energy electronic dispersion to the ex- 
perimental precision, like that of single-layer graphene 
0, 01- I n scanning tunneling microscopy/spectroscopy 
(STM/STS) measurements the Landau level quantiza- 
tion of the material in a magnetic field was found to be 
essentially that of single-layer graphene Theoreti- 
cally it has been shown that this approximate decoupling 
of different layers is due to a relative twist of the layers 
with respect to each other A rcnormalization of 

the electron velocity 0, Ell , van Hove singularities [3] , 
and interlayer transport [14j have been discussed as resid- 
ual effects of the interlayer coupling. 

In a recent STM measurement on multilayer epitax- 
ial graphene [l5[ a spatially modulated splitting A < 
lOmeV of the zeroth Landau level (LLo) was observed. 
In view of the above this finding is intriguing, since the 
states forming LLo of an isolated layer of graphene with- 
out electron-electron interactions are degenerate. There- 
fore, either the observed splitting of LLo is due to 
electron-electron interactions, or the interlayer coupling 
manifests itself prominently in this measurement. In 
many ways the experimental data favors the latter sce- 
nario. One such indication is the observation of a sub- 
lattice polarization of the split LLo: there are regions 
of space where the branch of LLo that has positive en- 
ergy A/2 appears to consist of wavefunctions localized 
on the A-sublattice of the graphene lattice, while the 
lower branch, at negative energy —A/2, is localized on 
the B-sublattice. The implied local sublattice-symmetry 
breaking has a natural explanation in terms of the in- 
terlayer coupling: the coupling to lower graphene layers 
generically induces a difference between the local envi- 
ronments of the two sublatticcs of the top layer in the 
material, which is probed in STM. This is illustrated in 
Fig. [1] for a stack of two graphene layers with a relative 
twist. There are regions where atoms on the A-sublatticc 
of the top layer are closer to atoms in the bottom layer 
than those on the B-sublattice of the top layer and re- 
gions with the reverse situation. A second conspicuous 



feature of the STM data is a spatial modulation of the 
splitting of LLo: the regions where LLo is split appear 
to be arranged on a hexagonal supcrlattice with a lattice 
constant I w 70 nm. It thus shares the symmetries of 
the moire pattern characteristic of the twisted graphene 
bilayer shown in Fig. [T] — another strong indication that 
the observed splitting is due to the interlayer coupling. 




FIG. 1: (color online) Moire pattern created by two graphene 
lattices with a relative twist. Top layer A/B sublattice 
atoms are shown as small blue/cyan (dark/light) spheres 
and connectors; bottom layer A/B atoms are shown as large 
red/yellow (dark/light) spheres. A region of AA alignment 
lies at the center, where each top-layer atom has a neighbor 
in the bottom layer. The AA region is surrounded by three 
AB- and three BA-aligned regions where atoms on only one 
top-layer sublattice have direct neighbors in the bottom layer. 
As a consequence, the sublattice-symmetry is broken locally. 

Earlier theory of the interlayer coupling in graphene 
multilayers did not predict the observed splitting of LLo ■ 
In Ref. [lit we therefore proposed a phenomenological 
theory, modeling the different local environments of the 
A- and the B-atoms of the top graphene layer by a "stag- 
gered" electric potential Vab that has opposite sign on 
the two sublattices. This model qualitatively accounts for 
the main features of the experimental data. In this Letter 
we present a microscopic theory of the interlayer coupling 
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in rotationally faulted graphene multilayers. We reduce 
the problem to an effective model of the top layer of the 
material, which is probed in many experiments, such as 
STM. In order to conveniently explain the rich spatial 
structure of the system illustrated in Fig. [TJ and observed 
in Ref. we formulate our theory in real space,as dis- 
tinct from prior momentum-space approaches [7, 11, 14 1. 
The resulting Hamiltonian reduces to the phenomeno- 
logical model of Ref. 15| in certain limits and it likewise 
reproduces the main qualitative features of the measure- 
ments. Our theory moreover allows us to test quanti- 
tatively whether the interlayer coupling can explain the 
experimental findings (l5j . The answer is affirmative: us- 
ing the commonly accepted tight-binding parameters of 
graphene multilayers our theory predicts both the mag- 
nitude of the observed splitting and its magnetic field 
dependence in very good agreement with experiment. 

We analyze the electron dynamics in a graphene layer 
"0" when coupled to a second layer "1," twisted by a 
relative angle 9 (9 = 0° for aligned honeycomb lattices, 
cf. Fig. [J), neglecting electron-electron interactions. The 
corresponding dynamics in multilayers at perturbati yely 
weak interlayer coupling, such as in the experiment [15|, 
are obtained by summation over all layers coupled to 
the top layer 0. Twisted graphene bilayers have been 
described before @, 0-12| by a local interlayer coupling 



Hamiltonian with parameters fitted to experiment [16j . 



dr* (0)t (r)r(r)* (1) (r) + h.c. (1) 



Here, the spinors collect the amplitudes for electrons 
on the two sublattices of layer j <E {0, 1}. The interlayer 
coupling r has contributions at wavevectors — 6*-% 
where b^ are reciprocal vectors of the graphene lattice in 
layer j [ijj] • The Fourier components of T quickly decay 
with increasing wavevector [ij 1(J 12 1. In this Letter we 
therefore neglect all but the zero wavevector component, 
setting r(r) =7, such that the distinction between com- 
mensurate and incommensurate interlayer rotations dis- 
appears. This approximation is valid for energies e 3> V, 
where V is set by the Fourier components of V that di- 
rectly connect K-points of the two layers [l2| . We take 
the limit < 6 -C 1, when V <C 7 (in the experiment 
[IH 9 w .25° and according to the estimate V ~ 9 2r y of 
Ref. [ijj this approximation is justified at all accessed 
energies). 

In our limit < 9 <C 1 a long-wavelength description 
is appropriate, where the isolated layers j are described 
by Dirac model Hamiltonians (we set % = 1) 

H {3) =v j drJ2^{r) [rr v ■ (-iV + eA(r))] ip^(r). 

(2) 

Here, er„ = (va x , a y ) is a vector of Pauli matrices, v = ± 
is the valley spin, — e the electron charge, and v the elec- 
tron velocity in graphene. We have included an external 



vector potential A to describe a perpendicular magnetic 
field B. Eq. ([2|) acts on the long- wavelength spinors ip^ u 
defined by *^(r) = Y, v u$ v (r)i>$ v (r). We write the 

Bloch functions u£l(r) - {£ p eacp[iK^-(r-T^)]}/y/3 
in the "first star approximation" appropriate for the in- 
terlayer coupling problem [13]. Here, p sums over the 
three equivalent Brillouin zone corners K^J that form 

the Dirac point of valley v [lj|] and t$} gives the position 
of an atom on sublattice fi S {A, B} within the unit cell 
of layer j. In the long- wavelength theory (which neglects 
inter-valley processes) the interlayer coupling reads 



ffi„t = / *^tf»t(r)i„(r)tf»(r) + /i.c. 



(3) 



with a matrix t whose long-wavelength components have 
wavevectors 5 K pv = (R$ — 1)K^J. Here, Re is a rotation 
around the z-axis by angle 9. Retaining only those long- 
wavelength parts of t we find 



(4) 



where terms of order 9 are neglected, while terms of order 
9Kr are kept as they may grow large. 

We next integrate out layer j = 1 in order to arrive 
at an effective Hamiltonian Hq S (uj) = Ho + 5Hq S (ijj) for 
the top layer j = 0, with 



SHf(uj) = H int (uj + V - H,)- 1 ^ 



(5) 



We include an interlayer bias V that accounts for differ- 



ent doping levels of the two layers [24j . In general, Hq is 



nonlocal in space and it depends on the energy uj. In the 
limit of a large interlayer bias, however, | V| 3> uj, 7, 9v/a, 
the sum uj + V — Hi becomes momentum- and energy- 
independent to a good approximation. The spatial non- 
locality and the energy-dependence of Hg S then may be 
neglected and Hq becomes a conventional Dirac Hamil- 
tonian ([2]) with a matrix potential 



6Hf = j rfr^tf)t(r) 



which we parametrize as 



t„(r)4(r) /(0) 



V 



$W(r), ( 6 ) 



Ur)4(r) 
V 



V oS {r) + vvea, ■ A off (r) + m off (r) v 2 a z . (7) 



The interlayer coupling in this limit generates effective 
scalar and vector potentials V eS and A GfT , respectively, 
and a mass term cx a z m eS v 2 that implies an effec- 
tive staggered potential V%g = m v 2 in locally Bernal 
stacked regions. It follows from Eq. (|4|) that SHq^ oscil- 
lates in space with wavevectors k = (Rg — 1)6, where b is 
in the "first star" of reciprocal lattice vectors of graphene. 
We plot 5Hq„ in the parameterization of Eq. ([7]) in Fig. 

m 
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Now turning to the experiment [15| we note that at 
large interlayer bias V our theory takes the form of the 
phenomenological Hamiltonian proposed in Ref. flij . It 
then intuitively explains the main qualitative features of 
the experiment: perturbatively in 7, the energy shift of 
a LLq wavcfunction ipo v in valley v is given by 



feo,^ = (^|*iC(w = 0)l^ 



(8) 



The unperturbed LLo wavefunctions are localized on in- 
dividual sublattices. Therefore, if SHq S included a con- 
stant staggered potential Vab > 0, with potentials Vab 
and —Vab for atoms on the A- and B-sublattice, respec- 
tively, a splitting A = 8eQ tV —\ — feo,i>=-i = — 2Vab be- 
tween sublatticc-polarized LLo states would result, as ob- 
served experimentally: Vab would increase the energy of 
the states ^0,1/=- 1 localized on the A-sublattice and de- 
crease the energy of the v = 1 states, localized on the 
B-sublattice. For the space-dependent V^ B = m cS v 2 of 
Fig. [2] that splitting is still present locally, around the 
extrema of m , at sufficiently large magnetic fields B, 
when the LLo wavefunctions fit well into the regions with 
extremal m off . Comparison of Fig. [2]with Fig. 5a of Ref. 
[TBI shows that the thus predicted spatial symmetries of 



oc \m 



eft 



agree with experiment. For large B the split- 
ting approaches lims-j-oo A = — 2VJ%. With decreasing 
B, as the wavefunctions become more extended, A gets 
averaged over maxima and minima of m eff and it is sup- 
pressed, also in accordance with experiment. 

The experiment of Ref. [l5| . however, was not done 
in the high bias limit. The fact that in the measurement 
[IH tunneling into LLo occurred only at a finite bias volt- 
age Vstm ~ 60meV between STM-tip and sample does 
indicate a doping of the graphene layers at the surface. 
The difference between the chemical potentials of the top 
layer and the layers below after screening is expected to 
be \V\ < Vstm ~ 60meV. However, the large applied 
magnetic field 4 T < B < 8 T corresponds to a large cy- 
clotron frequency uj c = \/2v/Ib [13]) where Ib = 1/VeB 
is the magnetic length: ui c rj 105 meV at B = 8T. In 
this experiment therefore |V| < uj c and H{f is not local 
on the scale Ib on which the wavefunctions vary. 

The experiment also indicates that it is the coupling 
between the top layer and its next-to-nearest layer (that 
is the third layer from the top) that produces the ob- 
served splitting. One concludes this from the observation 
that the dominant moire of the STM topography, most 
likely due to the coupling of the top layer to its nearest 
neighbor, has a much smaller lattice constant I pa 4nm 
than the supcrlattice associated with the splitting of LLo 
with I pa 70 nm. The estimates of the next-to-nearest 
layer coupling in the literature vary [li, T^- 20|, but there 
is a consensus that the coupling constant is 7 < 40 meV. 
The physics at the energies u) = ±A/2 pa ±5meV where 
the splitting of LLo occurs is thus described by Hq S at 
|w|,7 <C \V\ <C lo c - In this limit the effects of the inter- 
layer coupling are perturbative, which allows us to deal 





(c) Af 



(d) Af 



FIG. 2: (a) Effective potential V cS , (b) effective mass m off , 
(c) Af, and (d) Af of Eq. as functions of rO/a in grey- 
scale. Scale bars span a unity increment in r8/a. Note the 
expected sixfold and threefold symmetries of V eS and m eS , 
respectively. A cfI transforms as a vector under rotations. 



with the non-locality of -ffo analytically. We evaluate 
Eq. ||HJ| at |w|,7 <C |V| oj c in the appendix. In accor- 
dance with the intuition gained from the limit V — > 00 
of the previous paragraph, the resulting A is extremal 
in locally Bernal stacked regions and the wavefunctions 
are sublattice-polarized. The qualitative agreement with 
experiment thus carries over to the non-local theory. 

Now comparing our theory also quantitatively with the 
experiment we first take the limit of a large magnetic 
field, when the wavefunctions fit well into the Bernal 
stacked regions. The maximal splitting A max , reached at 
B —> 00 in AB- or BA-stackcd regions, can be extracted 
from Eq. ([20|) of the appendix by taking the limit 6^-0 
at fixed B. We find 



\V\ 



(9) 



in our approximations. Estimating 7 by 7 = 75 ss 
38meV given in Ref. [3] we find that |A max | rj 5meV 
for V ~ 40 meV. Considering the uncertainties in our 
knowledge of 7 and V, this agrees well with the experi- 
mentally observed value A sa 1 meV. 

We next quantify the magnetic field dependence of A 
in the regions with maximal A at B — > 00 (that is AB- 
or BA-stackcd regions) by expanding Eq. (f2"0)) asymptot- 
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ically for SKl B > 1: 



confined states 



-b,/b 




(10) 



2vSK 



V 



■ cos 



bB 



The crossover field £>» = 58K 2 /4e, where the exponent 
in Eq. (fTU)) becomes of order 1 and A starts to be ex- 
ponentially suppressed, evaluates to B» = 4.8 T for the 
interlayer rotation angle 8 = 0.25° of the moire pat- 
tern in the experiment of Ref. [l5| . Also that crossover 
field compares favorably with the experiment, where the 
splitting A disappears between B w 4T and B « 6T. 
Clearly therefore, the interlayer coupling can account for 
the main features of the splitting of LL reported in Ref. 
[lit a ^ so 011 a quantitative level. 

We finally discuss the influence of the graphene layers 
in the experimental sample that we have ignored so far. 
The coupling of the top layer to layers further away than 
the third layer from the top is negligibly small. The cou- 
pling 7 to the second layer, however, is not: 7 rj 0.4 eV 
18 1 . As mentioned before, the STM topography of Ref. 
15 1 has a moire pattern with scale I ~ 4nm, which in- 
dicates a rotation angle between the top two layers of 
8 m 4°. At this angle the coupling between the "first 
stars" of the Brillouin zones of those two layers is per- 
turbative, because of large energy denominators 0- The 
coupling between other X-points in the extended Bril- 
louin zone is too small to play a role at the scale of the 
observed splitting A [12]. The perturbative calculation 
outlined in the appendix therefore describes also the cou- 
pling between the top two layers of the measured sample. 
Applying Eq. p0[) to that coupling we find an exponential 
suppression of A that is lifted only above a crossover field 
B* = (9/8) 2 B^ « 3005* that is much larger than the ex- 
perimentally applied fields. The only interlayer coupling 
relevant to the experiment of Ref. [IH is therefore the 
next-to-nearest layer coupling discussed above. 

We conclude that the interlayer coupling is a viable 
explanation of the splitting of LLo reported in Ref. [l5| , 
both qualitatively and quantitatively. The theory that 
allowed us to reach these conclusions reduces in certain 
limits to an effective Dirac model for the top layer of a 
multilayer system, with effective potentials and a space- 
dependent mass. As such it makes the wealth of knowl- 
edge and intuition existing for the physics of single layer 
graphene available for the study of rotationally faulted 
multilayer graphene. Our theory thus appears to be an 
advantageous starting point for the exploration of much 
of the physics of this rather complex system. Numerous 
unconventional and so far unexplained phenomena ob- 
served in the material 2l| as well as known properties 
of our theory promise that such exploration will be re- 
warding. Especially the effective mass term is expected 
to have profound implications, for instance topologically 



Appendix: Perturbative Landau level splitting in a 
large magnetic field 

We evaluate the splitting A = Seo tV —i — <S£o,i>=-i be- 
tween the two valleys of LL at |w|,7 <C \V\ <§; w c , when 
it is perturbative, using Eq. (|8|) with localized wave- 
functions of LL : ipo, u =i = (0,cxp[— (a; 2 + y 2 )/4Z| + 
ixy /2l B ]/ 'V2ttIb) and ipo,v=-i = &yipo,v=i- We write the 
effective Hamiltonian as 



6Hf(r, r', to) = t v (r)G u (uj, r, r')tl(r% (11) 



where 



oj + V - se n 



n>0,s=± 



(12) 



with 



1 



s$„ (jL - kl B 



in terms of the oscillator wavefunctions 



(-1)* 



nl\/TT 



dx n 



e lkx (n>0) 
(13) 



(14) 



Here, $_i = 0, e n = \fnuj c and the wavefunctions in the 
valley v = — 1 are obtained as ip s ,n,k,v=-i = Oy^snkv=\- 
In our limit |V| <C w c , the contribution to feo.i/ with the 
smallest energy denominator comes from the term in Eq. 
(fT2|) with n = 0. That term is oc |tss| 2 in valley v = 1. In 
valley v = — 1 the corresponding term is identical, except 
that tBB is replaced by tAA- One has tsB = tAA + O(0) 
[25I ] . To leading order in 8 this term therefore does not 
contribute to A. The dominant contribution to fco thus 
comes from the off-diagonal elements of 5H cS and from 
the diagonal elements that are oc |^as| 2 or cx |£ba| 2 [the 
upper diagonal element in Eq. (| 12[) at v — 1]. In those 
matrix elements all contributing energy denominators are 
of the same order, O(io c ). We thus need to carry out the 
sum over n in Eq. (fl~2|) . We do this below for G\. The 
Green function in the other valley is then obtained as 
G-\ = a y G\(jy. We first rewrite Eqs. (fT2")) with (flU)) and 
dHl> as 

Gi(u,r,r') = f— e iHx-z')+i(y/iB-ki B ) 2 +( y '/i B -ki B ) 2 ]/2 



xg I u>, Mb, 7 Mb 

IB IB 



(15) 



and note that in our limit |a;| <C \V\ <C w c the component 
9AA-, which makes one of the leading contributions to A 
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according to the above considerations, can be expressed 
as 

9aa(0,x,X , ) = ^1(x,X , ) + o(J^\ (16) 

in terms of a function I that solves the differential equa- 
tion 

^ = -^i:*»W*»(x')e-<^"" 2 . (17) 

Using the completeness of the oscillator wavefunctions 
we find that Eq. (|17l) is solved by 

Kx,x') = [erf(x)-sgn(x-x')][ Grl (x)+ s g n (x-x')] 
+f(x) + f'(x') (18) 

I 



with arbitrary functions / and /'. Exploiting the sym- 
metries g(u,x',x) = 9(u,x,x') and g{u,-x,-x') = 

g(oj,X,x') that are implied by Eqs. (fl"2"|) and fp~5|) . one 
finds that f = f and that / is an odd function of x- Now 
noting that according to Eq. (TT^J) (3>n\G(oj)\$ m ) = for 
I! ^ m one concludes that / = [26j]. The off-diagonal 
matrix elements of g are found similarly. To leading order 
in V they read 

g AB (o,x,x') = ^Kx,x') + o^), 

g B A(0,X,x') = .9ab(0,x',x)- (19) 

Eqs. (jTTj) . (fT5l) . (fl6|) and (|T8|) allow us to evaluate feo,i^=i, 
Eq. ([5]), to leading order in 7, yielding 



<feo,*=i = 2^2 / 51 [erf (x) - sgn(x - x')] [erf (x) + sgn(x - x')] 

x {t BA (6K p )t BA (6K p ,) + (uj c /V2V) [t BB (6K p )t BA (SK p ,)(c ■ 5K p ,)* + t BA (6K p )t BB (SK p ,)(c ■ 5K P )}} 

x e -ll[{2{SK v , x -SK p , x ) 2 + (SK p , x +8K p , x ) 2 + (SK p , y -8K pl y ) 2 +2 l (SK p , :c +SK pl J(SK p , y -SK p , y )]/4:-l B [(c-5K p )x+(c-SK p ,yx'] ^ 

(20) 

I 



where scalar multiplication with c = (1, i) maps a vector 
a onto its counterpart c ■ a in the complex plane. Here, 
all wavevectors SK p are evaluated in valley v = 1. The 
energy shift <5eo.y=-i in the other valley is obtained as 
in Eq. (|20[) . but with t replaced by o y tay and SK P eval- 
uated in valley v = — 1. In Eq. ([9]) of the main text, 
that is in the limit of large B, only the first term in the 
curly brackets of Eq. ([20]) contributes and the resulting 
splitting A has extrema in regions where the layers are 



locally Bernal stacked and |£ba| 2 — Kab| 2 is extremal. 

In the limit B — > 0, when the wavefunctions become 
more and more extended and start averaging over several 
unit cells of the moire superlattice, the splitting of LLo 
decays to zero. In order to quantify this decay of A, Eq. 
(f2"Ul) may be expanded asymptotically in a large SKIb- 
Then Ssq^ is dominated by the terms with the weakest 
exponential decay in SKIb, which give 



, T/V ^ t BA (6K p )t BA (SK p/ ) + (io c /V2V) [t BB (6K p )t BA (SK p ,)(c ■ 5K p ,)* + t BA (8K p )t BB {8K p/ ){c ■ SK P )] 

x e -l 2 B l(SK p -SK p ,) 2 +SK p -SK p ,-i(5K p xSK p ,)-z]/2 ^1) 



at SKIb 3> 1. Here, z is the unit vector along the z-axis. 
Again all wavevectors SK p are evaluated in valley v = 1 
and (5eo,i'=-i in the other valley is obtained by replacing 
t with (Jytijy in Eq. (|21|) and evaluating 5K P in valley 
v = — 1. The sum over p and p' in Eq. (|21|) results in Eq. 
(fl"0| of the main text. 
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